Spontaneous creation of macroscopic flow and metachronal waves in an array of cilia 



Boris Guirao and Jean-Frangois Joanny 
Laboratoire Physico-chimie Curie (UMR 168), Institut Curie, 26 rue d'Ulm, 75248 Paris Cedex 05, Prance 

(Dated: February 2, 2008) 

Cells or bacteria carrying cilia on their surface show many striking features : alignment of cilia 
in an array, two-phase asymmetric beating for each cilium, coordination between cilia and existence 
of metachronal waves with a constant phase difference between two adjacent cilia. We give simple 
theoretical arguments based on hydrodynamic coupling and an internal mechanism of the cilium 
derived from the behavior of a collection of molecular motors, to account qualitatively for these 
cooperative features. Hydrodynamic interactions can lead to the alignment of an array of cilia. We 
study the effect of a transverse external flow and obtain a two-phase asymmetrical beating, faster 
along the flow and slower against the flow, proceeding around an average curved position. We show 
that an aligned array of cilia is able to spontaneously break the left-right symmetry and to create 
a global average flow. Metachronism arises as a local minimum of the beating threshold and leads 
to a rather constant flow. 



I. INTRODUCTION 

Many cells and bacteria have cilia or flagella on their 
surfaces. Examples are sperm cells which have one flagel- 
lum used for propulsion, the green alga Chlamydomonas 
that uses two flagella, and the much studied protozoan 
Paramecium which is covered by a layer of cilia. This 
layer is made out of approximately four thousands cilia 
which produce a very efficient motion with a velocity of 
order 1mm/ s in water, corresponding to 10 times the 
Paramecium size/s. Humans have ciliated cells in sev- 
eral organs: in the brain (cerebrospinal fluid flow), the 
retina (photoreceptor connective cilia), the respiratory 
tract (epithelial cells), the ear (hair bundles), the Falop- 
ian tube or the kidney... 

Cilia have two major roles: (i) detection (sensory cilia 
or flagella), for example in the retina, the ear and the 
kidney (ii) propulsion or creation of fluid flow (motile 
cilia or flagella) as for Paramecium or in the respiratory 
tract where the fluid flow is used to move away the mucus. 

The common structure of most cilia and flagella is an 
axoneme wrapped by the plasma membrane. The (9-1-2) 
axoneme is made of 9 microtubule doublets arranged on 
a circle around a central pair of microtubules 0. The 
cilium or flagellum is attached to the cell membrane by 
a basal body made of 9 microtubule triplets which has a 
structure very similar to that of a centriole. The basal 
body is attached to the cell membrane by anchoring fibers 
0- Typically the radius of an axoneme is 100 nm. The 
main structural difference between cilia and flagella is 
their length. The typical length of a cilium is lO/zm 
whereas a flagellum can be ten times longer. 

Dynein molecular motors are attached to the 9 mi- 
crotubule doublets; they move towards the microtubules 
— ends linked to the basal body and exert forces on 
the microtubules. Upon consumption of Adenosine- Tri- 
phosphate (ATP), dynein motion generates forces that 
induce a sliding between adjacent microtubules. Because 
the whole structure is attached at its basis, this sliding 
motion induces the bending of the cilium or flagellum 
and its beating. 



We here focus on ciliated cells creating fluid flow. 
These are cells with cilia on their surface, beating in one 
preferred direction in a coordinated way. One central fea- 
ture of cilia beating is the existence of two phases with 
a broken symmetry. Each beating can be decomposed 
into an effective stroke (ES) that propels the fluid and 
a recovery stroke (RS) where the cilium is coming back 
against the flow. In the example of Paramecium in water, 
the effective stroke lasts typically 9ms whereas the recov- 
ery stroke lasts 26ms. The typical beating frequency in 
water is 30Hz ,31]. The beating of Paramecium cilia is 
3-dimensional but for some species like Opalina the cilia 
remain in the same plane during their beating and the 
beating is 2-dimensional. In this work, we discuss the 
role of an external velocity field in this left-right symme- 
try breaking between the effective stroke and the recovery 
stroke for planar beating. 

One of the most striking features of an assembly of 
beating cilia is that they all beat in the same direction: 
the surrounding fluid can only be propelled efficiently if 
all the beatings have the same orientation. In all ma- 
ture ciliated cells, the beating direction is defined by 
the anchoring of the basal foot on the basal body. Only 
newly formed or developing cilia are randomly oriented 
01 . When they start beating, they tend to spontaneously 
align to finally beat in the same direction. One of the 
questions addressed in this article is the nature of the 
parameters that control this orientation. 

The role of the central pair of microtubules in the cen- 
ter of the axoneme is also a fundamental and complex 
question. In many species (such as Chlamydomonas), 
the central pair is both rotating and twisting within the 
axoneme during the axoneme movement. Current models 
postulate that the central pair modulates dynein activ- 
ity along outer microtubule doublets 5]. It thus allows 
the axoneme motion because if all the dyneins were act- 
ing at the same time, no bending would occur. Evidence 
in support of this model includes the observation that 
sliding between adjacent doublets occurs preferentially 
along doublets closest to one of the two microtubules of 
the central pair (the CI) in Chlamydomonas flagella j^. 
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Nevertheless, there exist also motile cilia with a (9+0) 
axoneme having no central pair. This means that cilia 
beating is possible even in the absence of the central pair 
of microtubules. Despite its importance, we do not dis- 
cuss the role of the central pair in the present work and 
we postpone its discussion to future work. 

Another important feature of ciliated cells, is the exis- 
tence of waves propagating all along the surface. These 
are called metachronal waves and might be due to the 
coordination of adjacent cilia for example via hydrody- 
namic interactions. Experimentally metachronal waves 
are observed to propagate in all possible directions: in the 
direction of the effective stroke (symplectic metachronal 
waves), in the opposite direction (antiplectic), or even in 
a perpendicular (laeoplectic or dexioplectic) or oblique 
direction. The origin of these waves and the mechanisms 
controlling their formation are not well understood. We 
show in this article that metachronism can arise naturally 
from the hydrodynamic couplings between cilia. Using 
a two-state model for the dynein motion as an internal 
mechanism of the cilia, metachronism appears to be a 
local minimum in the oscillation threshold of the motors 

m 

A last important feature of cilia beating that we wish 
to mention, is the role of calcium ions. The local [Ca^+] 
concentration has a strong influence on the beating pat- 
tern of cilia or flagella. For example detergent-treated 
Paramecium are able to swim forward at low [Ca^"''] con- 
centration (< 10~^M) and backward at high [Ca^+] con- 
centration (> lO^^M) because of ciliary reversal: the 
directions of effective and recovery strokes are switched 
[^Il0j|. In any case, the wild type Paramecium can have 
a very efficient backward motion monitored by calcium 
tanks in its body. We only discuss here qualitative as- 
pects of the role of calcium. 

In this paper, we address the question of the sponta- 
neous alignment of an array of beating cilia and the possi- 
bility of a spontaneous symmetry breaking in the beating 
that leads to the appearance of a macroscopic fluid flow. 
The internal mechanism of the cilia is described by the 
model of references 0, Hi which is based on a two-state 
model to describe the cooperative effects between dynein 
motors and only considers the relative sliding of two mi- 
crotubules in the axoneme. The coordination between 
the cilia is due to hydrodynamic interactions which are 
discussed in details in a coarse-grained description where 
the effect of the cilia on the flow is replaced by an effec- 
tive force. The outline of the paper is as follows. In the 
next section, we give a simple model for the alignment 
of beating cilia. In section IIIII we discuss the beating of 
one cilium following the model of Jiilicher and Camalet 
0, Hi ■ Finally, in section IIVI we discuss the spontaneous 
breaking of the left-right symmetry of the beating due to 
the flow created by the cilia themselves. 



II. SPONTANEOUS ALIGNMENT OF AN 
ARRAY OF CILIA: A SIMPLE MODEL 

A. Experimental results 

In an assembly of cilia covering the surface of a mature 
cell, cilia are beating in a preferred direction, and only 
newly formed or developing cilia are randomly oriented 
0. We first discuss the experiments showing how this 
preferred orientation is chosen. 

As mentioned before, the ciliary axoneme grows from 
a basal body analogous to a centriole. Two basal body 
appendages, the basal foot and the striated rootlet, lo- 
cated in the axial plane of the effective stroke, confer an 
asymmetrical organization to the basal body. The basal 
foot is laterally associated with two consecutive trip lets 
and points in the direction of the effective stroke [TllIT^ . 
The striated rootlet, associated with the proximal end of 
the basal body, sinks into the cytoplasm in the opposite 
direction . These two appendages define therefore an 
orientation of a cilium independent of the beating mo- 
tion. 

During ciliogenesis, newly formed basal bodies migrate 
toward the cell membrane where they anchor with no 
apparent order. Anchoring induces axoneme assembly, 
and cilia grow in random orientations. While cilia are 
growing, they do not beat immediately. A reorientation 
by rotation of the basal bodies in a common direction 
occurs at the final stage of ciliogenesis, when mature cilia 
beat ^3| ■ The preferred direction of the assembly is then 
well-defined. In the immotile-cilia syndrome, axonemes 
are incomplete, and the ciliary activity is abnormal or 
absent: the fluid is poorly or not propelled. On the cell 
level, the basal bodies are randomly oriented [T^ . 

These experimental facts suggest that the beating and 
orientation of cilia are closely related. Our working hy- 
pothesis is that the alignment of an assembly of beating 
cilia is mostly due to hydrodynamical coupling between 
cilia. The global flow created by the other cilia tend to 
orient a given cilium and above a certain beating ampli- 
tude, all cilia orient in the same direction. We now give 
a very simple modelling of this cooperative alignment. 

B. Alignment transition 

We assume in the following that the beating is planar. 
It is the case for Opalina for example, but not exactly 
for Paramecium where the recovery stroke is not in the 
plane of the effective stroke. 

Near the top of the ciliary layer, observations show 
that the velocity is time independent and uniform [l5l |. 
Consequently, we average the beating over one time pe- 
riod and replace each cilium of length L (and its effective 
and recovery stroke) by a single force (stokeslet) /, par- 
allel to the surface, created in the fluid of viscosity rj at 
height h < L above the membrane, as sketched in figure 

m 



3 




y 

/ 

-V- • — ► — ^ 










I 



FIG. 1: a) Beating pattern of a single cilium showing the 
Effective Stroke (ES), where the fluid is efficiently propelled, 
and the slower Recovery Stroke (RS), where the cilium comes 
back close to the surface to minimize the viscous effects, b) 
Effective force in the fluid / applied at a height h above the 
cell membrane, to mimic the cilium beating. 



We choose the x axis in the direction of the effective 
stroke and the z axis perpendicular to the cell surface 
that we approximate by a plane. The stokeslet along the 
X axis (/ = /e^), is located at point S — (0, 0, K). 

The velocity created at point X = {x,y^z) by this 
stokeslet with a no-slip boundary condition on the plane 
z = is given by : 

where the response tensor tensor G is given in |l6j | and 
reads: 



PiPk 

p3 
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RjRk 
RjRz 



(1) 



where p ~ X ^ S, R = X + S and a — x,y. 

In order to simplify the hydrodynamic problem, we 
assume in all the following that the cilia are far away from 
each other. In this asymptotic limit, it is consistent to 
describe the effect of the cilium in the fluid by a stokeslet. 
We introduce the two-dimensional vector along the cell 
surface r = (x, y) and consider the hmit r ^ z,h, L. We 
are interested in the velocity in the vicinity of the ciliary 
layer, typically 0<z<1.5L. At lowest order in z/r, 
the velocity reads: 

vir, d, z) = —er + 0{z^/r'') ~ u{r, 0) z (2) 

where we have defined a dimensionless height z = z/ L. In 
the following, we use mostly the velocity u(r, 9). At this 



FIG. 2: Square lattice of cilia with a distance d between two 
neighboring cilia. Cilium j exerts on the fluid a force / in the 
direction 4>j', Sji = (e^ifij 
cilium j to cilium i. 



order, the flow field is a radial flow centered on the cilium- 
stokeslet. Note that because of the no slip boundary 
condition on the surface, the force appears in the velocity 
field (Eq. [SJ in the combination fh homogeneous to a 
momentum. 

We consider now a regular array of cilia; this is a 
reasonable assumption for Paramecium, which shows a 
beautiful and very regular array of cilia on its surface 
We assume that this array is in an infinite plane. 
Cilium i is defined by its position in the x, y plane (the 
cell surface) by a vector t\ and by the angle of its plane 
of beating with the x axis, (pi as displayed on Fig. |2| 

The total velocity at the cilium i at height z, V{ri, z) 
is the sum of all the velocities Vj{ri,z) created by the 



other cilia j ^ i'. 



V{n,z) 



vj{n,z) 



We single out the the z dependence of V writing 
V{n, z) = U{ri) z with U{ri) = Y^j^i '^ji.^i) and 



Uj{fi) 



3fhL cos( 



2TT'ri 



(3) 



where eji is a unit vector from cilium j to cilium i and 
Oji = (el, eJi), as shown on Fig. El 

We now use a mean field approximation, replacing 
the velocity Uj{ri) by its average over the directions 
< Uj{fi) >0 given by: 



(4) 



where P{4') is the probability for a cilium to make an 
angle (jj with the x axis. The average velocity at cilium i 
is U{n) ~ J2J^^ < Uj{ri) >4,. 
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The mean field approximation assumes that the fluc- 
tuations around a given angle determining the direction 
of the flow are small. We choose the x axis in the direc- 
tion of the flow without any loss of generality, so that the 
probability P{(l>) is peaked around (p — 0. 

In order to determine the probability P{(f>), we write a 
stationary Fokker-Planck equation dtp J — 0. The prob- 
ability current J = Pdt<t> — Drd^pP is the sum of two 
terms: a convection term and a diffusion term where 
is a rotational diffusion coefRcient. The beating plane 
can fluctuate due to thermal fluctuations. Because of 
the flow, if the beating plane of one cilium is at an an- 
gle <f> with the flow direction, the cilium is subject to a 
torque M/'""" -aU sin0 along the z axis that tends 
to align it in the direction of the flow, where Uz is the 
velocity of the global flow and a a viscous coefficient in- 
volving the geometry of the cilium. A rotating cilium is 
also subject to a viscous torque jVf^^scous opposing the 
rotation M^^'^^"^'^ — —Cdt4> where ( is the rotational fric- 
tion constant. The total torque on the cilium vanishes 
(M/'"™ + = 0) and the probability distribution 

satisfies the Fokker-Planck equation 



d^P 



yU d 



[P sin (/)] = 



Defining the effective temperature as Dr( — ksT and 
imposing the normalization condition J P{(f))d(j) — 1, we 
obtain: 



P(0) 



(5) 



where /o(a;) is the modified Bessel function defined in 
[Tsj l . The average velocity U can then be self-consistently 
determined by calculating < Uj{fi) >0, and summing 
over all the lattice sites. We obtain 



(6) 



where K is a, constant depending on the nature of the 
lattice, d is the lattice constant (the distance between 
cilia) and Ii is a modified Bessel function For a 

square lattice, p^: 



K 



(-1,2 I l2\b/2 

(fc,;)#(o,o) ^ > 



2/3(^)C(^)^4.52 



(7) 

/3(s) and C(s) being respectively the Dirichlet and the 
Riemann functions 18] . 

The self-consistent equation for the flow velocity can 
be discussed by expanding the integrals Iq and Ii in the 
vicinity of [/ = 0: there are two solutions C/ = 0, and 
a solution at a finite velocity which exists only within a 
certain range of parameters 



U = 2\/2. 



ATrksT-qd^ 



SKafhL 



(8) 



This solution exists only if 

3KafhL 



AnkBTrid^ 



> 1 



(9) 



Within the mean field approximation Eq. |3 defines a 
dynamical phase transition between a non-moving fluid 
with randomly oriented cilia and a moving fluid with a 
global flow V{z) — Uz ^ given by Eq. |S1 where all 
cilia are spontaneously aligned in the same orientation. 
This dynamical phase transition is second order (with a 
continuous velocity at the transition) and it is associated 
to a spontaneous breaking of the initial 0(2) symmetry. 

The influence of some of the parameters can be di- 
rectly analyzed on Eq. A decrease of the distance d 
between two cilia favors the alignment, increasing the hy- 
drodynamic coupling. A decrease of the temperature T 
also favors alignment as the random thermal motion op- 
poses the alignment. An increase of the effective hydro- 
dynamic force of one cilium / is associated to an increase 
the hydrodynamic interactions between cilia and leads to 
a better alignment. The same effect occurs for a and the 
cilium length L. Finally, increasing h helps to create a 
global flow, since the velocity on the membrane vanishes 
and the higher the force is exerted, the more efficient. 

A more precise analysis requires the estimation of the 
parameters /, a and h. The height is of the order of the 
cilium size h ^ L. The calculation of a is given in ap- 
pendix I for a general beating (Eq. I54|l . It turns out that 
a is linked to the difference of the areas covered during 
the effective and recovery strokes. Here we approximate 
a ^ (^±LA, where is the perpendicular friction con- 
stant per unit length of the cilium and A the amplitude 
of the movement of the tip. 

In order to give a simple estimation of the effective 
hydrodynamic force, we consider that the friction coeffi- 
cient takes the perpendicular value during the effec- 
tive stroke, and the parallel value ^|| during the recovery 
stroke. Introducing the beating frequency uj, we estimate 
/ ^ i^-L ~ £,\\)LujA. a more precise calculation of these 
two quantities as a function of the beating patterns is 
given in appendix I. 

Consequently, we obtain: 

SKafhL (.±{£.± - e||)-A2Lcj 



AnkBTrid^ 



ksTr) 



d3 



The difference between the two local drag coefficients 
and is the key to an efficient beating. Increasing the 
amplitude A or the frequency u of the beating favors 
the alignment of cilia, as could be expected. Increasing 
the viscosity of the medium also promotes the transition, 
because it increases the coupling between cilia. 

This naive mean field approximation is only a first step 
of our study. In the following, we take a closer look at 
the internal beating mechanism of one cilium and then at 
the beating of an array of cilia to obtain a more precise 
and quantitative description. 
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III. AXONEMAL BEATING 

In this section we discuss the beating mechanism of a 
single cilium. We follow closely the work of Camalet and 
Jiilicher which mimics the cilium by two microtubule 
filaments sliding along one another under the action of 
the dynein motors and uses a 2-state model to describe 
the collective motion of the dyneins. We use as bound- 
ary conditions for the motion those introduced recently 
by Hilfinger and Jiilicher |23| that seem to have good ex- 
perimental support |2ll |. In the next section, we use the 
same model to discuss the coordination between cilia. 



A. Equation of motion 

Each microtubule doublet within the axoneme can be 
described effectively as an elastic rod. Deformations of 
this rod lead to local sliding displacements of neighbor- 
ing microtubules. Here, we only consider planar defor- 
mations. In this case the geometrical coupling between 
bending and sliding can be captured by considering two 
parallel elastic filaments (corresponding to two micro- 
tubule doublets) with a constant separation a along the 
whole length of the rod (see Fig. EJ. At one end, which 
corresponds to the basal end of an axoneme, the two fil- 
aments are elastically attached and are allowed to slide 
with respect to each other, but not to tilt ^J. The basal 
connection is characterized by an elasticity k and a fric- 
tional drag 7. The configurations of the axoneme are 
described by the shape of the filament pair given by the 
position of one filament X(s) at arclength s. The shape of 
the other filament is then given by X'{s) — X{s) — an{s), 
where n is the filament normal. In the following, we de- 
scribe the filament conformation by the angle ■0 between 
the local tangent vector and the z axis or by the defor- 
mation h in the transverse direction defined in Fig. 1111 

The energetics of the filament pair is due to the bend- 
ing elasticity. In addition to filament bending, we also 
take into account internal stresses due to the active ele- 
ments (dyneins). We characterize them by the force per 
unit length /(s) acting at position s in opposite directions 
on the two microtubules. This force density corresponds 
to a shear stress within the cilium which tends to slide 
the two filaments with respect to each other. 

The local curvature is C = dgtp (see Fig. Illll . The slid- 
ing displacement, A{s,t) is related to the the sfiding dis- 
placement at the base Ao(i) by A(s, t) — Ao{t) — atp{s, t), 
because we impose the boundary condition ip{0,t) — 0. 

A configuration of a filament pair of length L is asso- 
ciated to the free energy functional: 

Here, k denotes the total bending rigidity of the fila- 
ments. The inextensibility of the filaments is taken into 




FIG. 3: Two filaments (full curves) X and X' at constant 
separation a are rigidly connected at the bottom end where 
s = 0. Internal forces /(s) are exerted in opposite directions, 
tangential to the filaments. The displacement A at the tip is 
indicated. 

account by the Lagrange multiplier A(s) which enforces 
the constraint (dgX)'^ ~ 1. 

The first term of this equation is the elastic energy 
due to the basal sliding occurring with a connection of 
elasticity fc. 

The tangent component of the integrated forces acting 
on the filament between s and L is denoted by r(s). As- 
suming that there is no external force applied at the end 
s = L oi the cilium: 

r(.)=r(s). fds' ^-^ = -t{s). rds'^ 
Js oX Jo oX 

We assume for simplicity, that the hydrodynamic ef- 
fects of the surrounding fluid can be described by two 
local friction coefficients and ^|| for normal and tan- 
gential motion. The total friction force per unit length 
exerted by the cilium on the fluid is then fy[X{s)] — 
{£,\\t t + ^±n n) dtX{s). The force balance at arclength s 
can then be written as 

1 1 SG 

dtX = ~(—n n + —ti) — (10) 

C-L ^11 SX 

which leads to: 

dtX = —{-K'li) -af + i)T) + —{ntPij + aipf + i) (11) 
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where the derivatives with respect to arclength have been 
denoted by a dot. 

The beating of the filaments is very sensitive to the 
boundary conditions imposed at its ends. As the the 
force density —6G/SX is equiUbrated by the density of 
friction force exerted by the fluid on the system (see Eq. 
llOfl . the boundary contributions coming from the free 
energy variation SG arc equihbratcd by external forces 
fext and torques Text — Text (iy applied at the ends. 

At the free end of the cilium, both the external force 
and the external torque vanish and 



fext = -[nC + af'jn + Tt = Q, Text ^ ^ (12) 
At the base, s — 0, the boundary conditions are: 

fext = (kC* + a/j n - Tt, Text = -kC + a j dsf{s) 

(13) 

The external torque and force are chosen in such a way 
that the base is fixed {dtX = 0) and that the cilium 
remains perpendicular to the surface ( dtt = Q ov il>[Q) — 
0). 

The final boundary condition is associated to the basal 
sliding 

7atAo = - -fcAo + r dsfis) (14) 

The determination of the cilium motion requires a 
model for the shear force created by the dyncins that 
we calculate using a two-state model. 



B. 2-state model for the cilium 

Following 0, we now introduce the 2-states model of 
coupled molecular motors [23, 113 to describe the internal 
mechanism of the cilium. This model allows the calcula- 
tion of the shear force / due to the dyneins. 

Each motor has two different chemical states, a 
strongly bound state, 1, and a weakly bound state, 2. 
The interactions between a motor and a filament in 
both states are characterized by potential energy land- 
scapes Wi{x) and W2{x), where x denotes the posi- 
tion of a motor along the filament. The potentials have 
the filament symmetry: they are periodic with period /, 
Wi{x) = Wi{x -\- 1) and are, in general, spatially asym- 
metric, Wi{x) 7^ Wi{~x). 

In the presence of ATP, the motors undergo transitions 
between states with transition rates oji and uj2- Introduc- 
ing the relative position ^ of a motor with respect to the 
potential period, [x — + nl with < ^ < ^ and n an 
integer), we define the probability Pi{(_,t) for a motor 
to be in state i at position ^ at time t. The relevant 
Fokker-Planck equations are: 



dtPi +v d^Pi 
dtP2+vd^P2 



-UJlPl + UJ2P2 
UJlPl - UJ2P2 



where v = dtA = adtipis) is the sliding velocity between 
the 2 filaments. 

The simplest choice of the two potentials Wi, is a saw- 
tooth potential (with barrier height U 3> kT) represent- 
ing a strongly bound state for Wi , and a flat potential W2 
representing a weakly bound state. Here, for simplicity, 
we use the symmetric potentials: 

Wi{x) = J7sin2(7ry) 



W2{x) 



W2 



Although this choice is somehow arbitrary, we checked 
that the final results only depend qualitatively on the 
actual shape of the potentials and of the rates defined 
below. 

When a number of motors act together to propel a 
filament, however, the direction of motion is a collective 
property. The filament might move in either direction 
|23| . The absence of asymmetry in the potentials implies 
that an individual motor is not able to move directionally. 
It is not the case for an assembly of motors: even with a 
symmetric potential, provided that detachment can only 
take place at a localized position near the bottom of a 
potential well, oscillations can occur. 

We define the distance from equilibrium Q: 



it = Sup[oj[ I e 

UJ2 



I oc 6^'^/'=^ - 1 (15) 



is related to the chemical potential difference between 
ATP and its hydrolysis products, A/i — ^atp — (J'Adp — 
lip. At equilibrium, A/i = and = 0. We assume for 
simplicity that the binding rates W2 and the detachment 
rate Wi are given by: 

UJ2iO = '^(l + ^sin^TT^)) 



Note that, with this choice the sum uji + uj2 ^ ^^(1 + ^) 
does not depend on ^, and that if SI = 0, (jJi = and no 
directional movement is possible. Here v is a constant 
transition rate. 

If we assume that the motors are uniformly distributed 
along the filaments with a density p, the probabilities Pi 
and P2 satisfy the relationship Pi + P2 = p. The Fokker 
Planck equation reduces then to a single equation for 
P^Pi 



dtP + idtA)d^P =-iLUl+ LU2)P + PUJ2{0 



(16) 



This model leads to an expression for the shear force 
per unit of length /(s, t) created by the dyneins and driv- 
ing the the cilium beating. Using the results of |23| , and 
the fact that W2 is a constant: 



1 

/(s, t) = -KA - XdtA -jj^d^ PiOdiWi 



(17) 
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where K is an elastic stiffness per unit length mimicking 
the influence of the nexins which are proteins acting as 
springs in the axonemal structure, and A is an internal 
friction coefficient per unit length modelling the friction 
encountered by the motors. Eauations llll and llGI allow in 
principle a complete calculation of the beating motion. 

In the following, we assume that the beating oc- 
curs with a "small" amplitude which means that both 
tp{s,t) ^ 1 and h{s,t) <^ L. A quick look at the beating 
pattern of a cilium of Paramecium shows that the beat- 
ing occurs with a large amplitude. Nevertheless, this ap- 
proximation allows us to extract interesting information 
on the parameters controlling the beating. Moreover, 
the work of Hilfinger shows that larger amplitude beat- 
ing patterns are very similar to small amplitude patterns 
[20|. We must however keep in mind that our approach 
is valid only if the the system stays close to an oscilla- 
tion bifurcation which is consistent with the fact that we 
consider only small movements. 

We use the deformation h, rather than V' or A to de- 
scribe cilium motion and work at second order in \h\ so 
that ip = h + 0{\h\'^). In the absence of any external 
flow, the equation of motion 1111 projected on t imposes 
that r = 0(|/ip). The projection of the equation of mo- 
tion on n then yields: 



^^dth = -K'h -af + 0{\h\^) 



(18) 



The non-linear terms are not important here but they 
will turn out important in the following section. Indeed, 
experiments show that the beating is clearly asymmetri- 
cal (ES and RS), and we must expand at least to next 
order if we want to capture this phenomenon. With this 
variable and at this order, the boundary conditions read: 

/i(s = 0) = 0, A(s==0) = 
Kh{s = L)+ af{s = L) = 0, h{s = L) = (19) 

The explicit solution of the equation of motion ITTjl 
is obtained by Fourier expansion in time h(s, t) = 
^nL-oo ^n(s)e™"*. The explicit derivation of the equa- 
tion satisfied by the Fourier components is given in ap- 
pendix II. At linear order the effect of the motors is char- 
acterized by a susceptibility 



x{^, i-^) = — Xiuj + 



(20) 

Using dimensionless variables, s = s/L, u = oj, 

K 

Xn — x(f^j ™-^) and h = h/L, the equation of motion 

K 

for the Fourier component n reads 

h „ + Xn hn + inCu /i„ = (21) 
The boundary conditions at the base s = are 

hn{0) = 0, flniO) = (22) 



At the free end of the cilium, s = 1 : 

'Kn{l)+Xnhn{l)+fnhn{l)^0, ^„(1)=0 (23) 

with r„ = Xn/ (fc ~ Xn + in^iij) where we have intro- 
duccd the dimensionless parameters k = k and 

K 



1 ■ 
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C. Beating pattern 

In the absence of any external flow the beating is sym- 
metric and the Fourier components n=0 and n=2 of h 
vanish for symmetry reasons. Close to the oscillation bi- 
furcation threshold, cilium beating is dominated by the 
first Fourier component n = 1. 

The solution of the linear equation of motion |^ is 
written as the sum of 4 exponentials 

hi{s) = Ai (e«i"^ + bie-^^^' + cie'i^' + die''^^') (24) 

where the two inverse decay lengths are given by 



91 
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1/2 



1/2 



(25) 



The boundary conditions are explicitly discussed in ap- 
pendix II. The condition for existence of non vanishing 
solutions is given by equation 1631 of Appendix II. This is 
a complex equation, that gives therefore two conditions 
which determine both the the critical value of the dis- 
tance from equilibrium where the oscillations start flc, 
and the reduced oscillation frequency ujc- The critical 
value flc is a Hopf bifurcation threshold: there are no 
oscillations if 17 < fie and cilium beating is only possible 
if > Oc. At the bifurcation threshold, the amplitude 
of the oscillations vanishes. It is not possible to calculate 
the amplitude of the oscillations above the bifurcation 
threshold with the linear theory presented here. This re- 
quires a complete determination of the third order terms 
in the equation of motion which goes far beyond the scope 
of this work. This very complex problem is attacked in 
the work of Hilfinger and Jiilicher . 

An analytical determination of the Hopf bifurcation 
threshold and of the beating frequency at the threshold 
does not seem to be possible analytically. We therefore 
rely on a numerical solution, choosing reasonable values 
of the various parameters. 

We study a cilium of length L = 12/im, which is the 
length of a Paramecium cilium. It moves in a fluid of 
viscosity 77 ~ ^rj^ater — 4 \0~^Pa.s, which is higher than 
the water viscosity in order to take into account the pro- 
teins above the cell body. We estimate k = 4 XQ^^^J.m 
corresponding to 20 microtubules The other pa- 

rameters are I = lOnm, a = 20nm, U = lOkT, K ~ 0, 
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FIG. 4: Approximate cilium deformation h{s,t) at different 
times steps (corresponding to different colors) during a beat- 
ing period. Tlie beating is symmetrical with respect to the 
vertical axis. Deformations are propagating from base to tip. 
With ^1 = 1/70 the maximum deformation is hmax — 0.14. 



= 5 lO^TO^^, and A = IPa.s very similar to those of 
[7]. In order to match the typical frequency observed 
in Paramecium, and to obtain a realistic pattern of the 
beating, we take = 35.577 = 142 lO^^Pa.s and we 
choose k = 6.547Vm^, 7 = 7.31 -qwater, and v = 600s~^. 

The value of is rather high, but it must include 
the hydrodynamic interactions of the cilium with the cell 
surface (see ^2^) which are not taken into account if we 
use the classical friction per unit length of a rod. 

The numerical resolution of Eq. then yields 

(ric,<^c) — (6.55 10~^, 1294). This corresponds to a crit- 
ical beating frequency fn ~ 28Hz which is the typical 
value for Paramecium [Sj. The calculated beating pat- 
tern of the cilium is shown on figure 01 It corresponds to 
a wave (or a superposition of waves) propagating from 
the base of the cilium to the free tip as observed experi- 
mentally. 

A detailed study shows that the equation giving the 
bifurcation threshold and the beating frequency has sev- 
eral solutions. (ni"\w^"^) with > A first 
guess would be that the axoneme starts beating at the 
lowest threshold uii^"^ ) . However, it is known exper- 
imentally that during the beating the deformation waves 
propagate from the base to the tip and not from the tip to 
the base |2^ . The first two oscillating modes correspond 
to waves propagating in the opposite direction. In order 
to be consistent with the experiments we do not consider 
them here. A better choice of the transition rates lui and 
UJ2, would perhaps allow to justify this choice. The di- 
rection of propagation of the wave is extremely sensitive 
to the boundary conditions. We have allowed here basal 
sliding as suggested by some experiments and we have 
imposed that the cilium is clamped at its base with an 
angle ip = 0. This also seems consistent with some ex- 
periments analyzed in j20j |. The other extreme limit of 
a completely free cilium (a vanishing external torque at 
the base) leads to a wave propagating form the base to 



External viscosity Critical frequency fc Simulations 





28Hz 


29Hz 




WHz 


ITHz 




UHz 


UHz 



TABLE I: Decrease of the beating frequency with increasing 
external viscosity as observed in experiments |27l| . Compari- 
son with the simulations done in for one single cilium. 



the tip for the first mode. We have tried to use an in- 
termediate boundary condition where the torque at the 
base is an elastic torque and varying the related stiffness, 
however we were not able to obtain a beating pattern 
looking like the experimental one. We therefore proceed, 
considering only the third beating mode. 

The beating frequency fc varies with the viscosity of 
the medium. Experimentally, when methyl-cellulose is 
added in water, the viscosity increases significantly. We 
predict here a decrease of fc with increasing external vis- 
cosity as observed in the experiments of |27j and in nu- 
merical simulations (see Tabled). We observe an ap- 
proximate linear decrease of the beating frequency when 
plotted against log{ri/r]u]) as in the simulations performed 
in |23| (we find similar values for the frequency). 

The effect of [C'a^"''] on the beating pattern can also 
be studied qualitatively. As mentioned in the introduc- 
tion, [Ca^+] has a strong influence on the beating pat- 
tern. Calcium concentration variations are at the ba- 
sis of the shock responses of many organisms, changing 
the ciliary-type beating into a flagellar-type beating in 
Chlamydomonas, or switching the directions of effective 
stroke and recovery stroke in Paramecium, or in reversing 
the direction of the " wave" propagation on the flagellum 
and thus reversing the direction of the movement in Chri- 
tidia "2^. 

This last example can be explained qualitatively within 
our approach. In Chritidia, both directions are possible 
for the deformation wave propagation. Calcium may af- 
fect the chosen mode of beating, allowing the system to 
choose Vl'^P and its tip-to-base pattern instead of Jli^"* . 

On the other hand, calcium is also likely to change the 
attachment/ detachment rate (and thus change the pa- 
rameter v) or the boundary conditions at the base of the 
cilium (and thus change k). In Chlamydomonas, calcium 
has a contractile effect on the striated fibers connecting 
the basal bodies of the two fiagella 30] . These changes 
induce a change in the beating pattern, and may result 
in a switch from base-to-tip to tip-to-base wavelike prop- 
agation. 



IV. LEFT-RIGHT BEATING SYMMETRY 
BREAKING 

In the presence of a transverse external flow, the beat- 
ing can no longer be symmetrical as sketched on Fig. 
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FIG. 5: Effect of an external flow V on the beating of a single 
cilium. a) Symmetrical beating, b) broken symmetry due to 
the external flow. 



13 The cilium tends to beat faster and quite straight in 
the direction of the flow, whereas it comes back slower 
and more curved against the flow. This looks hke a two- 
phases beating with an effective and a recovery stroke. 

If the beating is asymmetrical, the cihum exerts a force 
in the fluid that can itself produce a flow. In a cer- 
tain range of parameters, one can therefore expect that a 
continuous flow is spontaneously generated by hydrody- 
namic interactions between cilia: an assembly of ciha, 
beating symmetrically, is able to break spontaneously 
this left-right symmetry of the beating to create a global 
flow. This idea of a spontaneous breaking spontaneously 
of the left-right symmetry has already been suggested in 
|3l| with a more abstract system (called rowers) having 
two internal energy states. 

In this section, we first study the effect of an external 
velocity imposed by the experimentalist on the beating 
symmetry of a single cilium. We then consider an ar- 
ray of aligned cilia and determine the conditions under 
which this assembly of cilia breaks its left-right symme- 
try and generates a global flow. Metachronal coordina- 
tion between cilia naturally emerges from hydrodynamic 
couplings as a local minimum of the oscillation threshold 



The equation of motion [TTl reads then: 

-* 71 •■■ •■ t ' " 

dtX = y-F — (-KV-a/+Vr)-K — (K7/;V+aV'/+T) (27) 

The boundary conditions are the same as in the absence 
of the neighboring cilia and are given by Eq. ^1 El 1141 
Following the same procedure as for a cilium in the 
absence of flow, we find the equation of motion for the 
deformation of the cilium h: 



Udth = - n'h -af- ^\iV{hh 



2^11 



+ Oi\hf) 
(28) 

The introduction of the external fiow breaks the h — > 
— h symmetry (or left-right symmetry) introducing in ll8l 
terms of zeroth and second order in h in the equation of 
motion. The boundary conditions do not depend on the 
external fiow. 

As above, we expand the deformation of the cilium h in 
Fourier components in time. Using the same notations as 
before, the equation of motion of the Fourier components 
can be written as 

_ ^ _ _ _ 17 _ ^2 

hn+Xn hn + inuj hn = F(5o,„ - —{^ hh + h )„ (29) 

for n = 0,1,2 and where we have introduced the new 
dimensionless parameters: 



In the limit of small external velocities, we have neglected 
terms of order . 

The equation for the first mode is identical to Eq. 1611 
with the same boundary conditions. At this order in 
V, the fundamental mode is not affected by the exter- 
nal fiow. Consequently, the oscillation threshold and the 
beating frequency are the same as in the absence of fiow 
and the Fourier component hi is given by Eq. 1241 

The zeroth Fourier component ho gives the average 
deformation of the cilium. It is a solution of 



External breaking of the beating symmetry: 
cilium submitted to an external flow 



We impose an external fiow V = Ve^ along the x axis 
for simplicity. It is found experimentally that the velocity 
above the cilia sub-layer is time independent and uniform 
|15|. justifying our choice. This fiow is in this first part 
externally fixed and we consider the limit of vanishingly 
small flows. 

The force per unit length exerted by the cilium on the 
fluid fv[X{s)] depends on the external velocity V. 



fAX{s)] = i^iitt + Un H) idtX{s) ~ V) 



(26) 



K hn^V 



V 
"2 



Uhihl + hlh) + \h\ 



(30) 



with the same boundary condition as before. Neverthe- 
less, ho does not vanish at first order in velocity because 
of the broken symmetry due to the external fiow which is 
refiected in the right hand side of Eq. 12111 The complete 
solution for ho is rather tedious to obtain and lengthy. 
We do not display it here explicitly. We write it as the 
sum of two contributions; h^'^, corresponds to the cur- 
vature of the cilium under the flow V at equilibrium, i.e. 
in the absence beating (hi = 0), and h-o^^ , corresponds 
to the corrections to this equilibrium deformation due to 
the beating when there is enough ATP in the medium 
ft,Q If as above, we ignore the elasticity of 
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a) 



b) 



FIG. 6: a) Average position of a cilium which is curved 
in the direction of the flow, ho(s) =< h{s,t) >. b) Second 
Fourier component of the deformation 25R[/!,2 (s)e^'"^*] at dif- 
ferent times during a beating period. The scale is dilated: 
\h2\ <C 0.1 with the parameters ^ = 1, Ai = 1/70 and V = 1. 



the nexins {K — > 0) 



V 



-ATP 



is) 



V 



(31) 



where Ai is the amplitude of the first Fourier mode of 
the oscillation defined in Eq. 1241 and 0o(s) is a linear 
combination of exponentials. In the limit V = 0, ho = 
as expected. The average deformation of the cilium is 
plotted on Fig. which shows the bent shape under the 
action of the external flow. 

The second Fourier component gives the asymmetry of 
the beating. It is obtained from the equation of motion 



^ 2 + X2 h2 + 2iuj h2 



V 
"2 



£_ hihi + hi 



(32) 



We do not give here the lengthy explicit expression of 
h2 but we write it as 

Ms) = jAlMs) 

where 02 (s) is a linear combination of exponentials. Here 
also, in the limit V = 0, h2 = 0- The plot against s at 
different times on Fig. ^ leads to a complicated pattern. 

The total deformation of the cilium h(s,t) ~ ho{s) + 
/ii(s)e'"'=* + /i2(s)e^*"'^* + c.c. is plotted against s at dif- 
ferent times equally spaced on Fig. [7| In order to stress 
the fact that the beating is easier and faster in the direc- 
tion of the flow, and more difficult and slower against the 
fiow, we have chosen rather large values of the parame- 
ters, ^ = 1, Ai = 1/5 and V = 2, and we plot h{s,t) for 
s G [0, 0.2] on Fig. 13 

The external fiow thus breaks the left-right symmetry 
in two ways. First the average position of the cilium is 



FIG. 7: Beating pattern at the basis of the cilium (s £ [0, 0.2]) 
with the parameters ^ = 1, ^i = l/5 and V — 2 : the cilium 
beats faster in the direction of the flow and slower in the 
opposite direction around a curved average position. 



not the vertical axis but a cilium curved in the direction 
of the fiow. Second, the beating itself is no longer left- 
right symmetric: the cilium goes faster in the direction 
of the flow and comes back slower against the fiow. The 
beating pattern looks like a two-phases beating with an 
effective stroke and a recovery stroke. The external flow 
may therefore be an important factor in the asymmetry 
of the beating. 

Another important result, is that, because the beating 
propagates a base-to-tip deformation, the curved cilium 
exerts a finite average force in the fiuid in the direction 
of the flow. Thus, if an external flow breaks the left-right 
beating symmetry, the cilia create a force in its direction 
and can amplify this fiow. This is the basis of the left- 
right spontaneous symmetry breaking that we discuss in 
the next section. 

The external flow is not always the only source of sym- 
metry breaking. If it were so, a Paramecium would al- 
ways go in the same direction once it started moving. 
This is not the case, this organism is able to go backward 
when it bumps into an obstacle thanks to the release of 
calcium that reverses the beating. 

The calculations of this section have been made with 
a velocity V uniform over the cilium length. This is 
not consistent with the presence of a cell wall where the 
cilium is anchored. Nevertheless, the main idea was to 
study how an external flow can break the beating sym- 
metry in the simplest way. Similar calculations can be 
performed with a linearly varying velocity V = Uz; they 
do not lead to any new physical effects. 



B. Spontaneous breaking of the beating symmetry: 
array of aligned cilia 

We now consider a regular array of cilia on a cell body, 
beating all in the same direction. Starting from a sym- 
metrical beating, we show that the left-right symmetry 
is spontaneously broken within a certain range of the 
parameters controlling the beating due to the hydrody- 
namic couplings between cilia. 
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1. Equations of motion 

For a cilium located in the xy plane at position r, we 
call y[X(s)] the velocity created by the other cilia at the 
point X{s) of arclength s. The equation of motion of 
the cilium is similar to that obtained previously with an 
external flow field and we write up to third order in h as 



The Fourier components of the force fo =< fx > and 
/i are calculated using the expression of and its 

average over one time period given by Eq. EH in the 
small movements approximation: 



fo ~ 2lj{^^ - ^\i)^[2hohihl ^ ill 
fi ~ iuS,±hi 



duhi{u)ho{u)] 



(37) 



^j^dth = -Kh -af + ^^n.V + 0{h^,t.Vh) (33) 

where the projection of the local external velocity on the 
cilium normal is 

rl.V = V^{l-h^/2-Vzh + 0{h'^)) (34) 

The boundary conditions for the motion are the same as 
in the previous section. 

The velocity Vj[Xi{si)] created at arclength Si of the 
cilium i by a cilium j is given by 



Vj[XM)]^ / ds, G[X,is,),X,{s,)].f,[X,{sj)] 
Jo 

where G is the second order hydrodynamic tensor given 
by Eq. and fj — f^'^"-* is the force per unit of length 
created by the beating of the cilium j. The total velocity 
at the arclength Si of the cilium i is thus given by 

V[XM)] ^J2vj[X,{s,)] = V{n,s,,t) (35) 

As in section m we consider the limit L <^ d and we 
only keep terms of the second order in s/r, r being the 
distance between cilia so that 

G./ = ([G./J.e;) e; + 0{s^ /r^) 

This means that only the velocity along the x axis created 
by the component fx of / plays a role and that we can 
ignore the other component Vz of the velocity. Using the 
notations of Fig. |51 and noting that z = s -\- 0{h^), we 
obtain 



where 5 is the imaginary part of a complex number. 

We assume that all cilia are identical, and that they 
all beat with the same pattern. The only difference in 
the beating patterns of cilia j and i is a possible phase 
difference that we call tpij . Defining 



and dropping the index i, we write the Fourier compo- 
nents of the velocity as 



Vi{s) 



SIqs ^ cos^ 9ij 



27r77 ^ \n - r 

3/is cos^ 9ij 
27r?7 ^ \n-rj\^ ^ 



3/C/o 



(38) 



3IC[{ifiij}]Ii 
2Trrjd'^ 



Vxin, Si,t) 



3si \ ^ cos^ dij 



^,3 / ds,f,^is,)s,+Oih\-) 
r-ii Jo ^ 



The geometrical constant /C is given by Eq. \7\ioT a square 
lattice of cilia spaced by d. The constant /C[{(/9y }] de- 
pends on the relative phases between the cilia. If the 
phases ipij are randomly distributed, then /C[{(^y }] ~ 
and Vi = 0. There is no oscillating component of the 
velocity. 

On the contrary, because we know that metachronism 
occurs in an array of beating cilia, we choose a constant 
phase difference ip between two consecutive cilia in the di- 
rection of the plane of beating : (fij+i — = f- This is 
the case for simplectic and antiplectic metachronal coor- 
dination. We only consider those cases (and not laeoplec- 
tic or dexioplectic metachronism) here. Experimentally, 
for Opalina (simplectic) and Pleurobrachia (antiplectic) 
that both have planar beatings, no metachronal wave in 
the transverse direction of the beating can be seen . 

We stress that we do not impose the phase difference 
If. The system is free to adjust its phase. We then write 
/C[{(^ij}] — IC{ip) with 



(36) 

As in the previous sections, we expand the velocity, 
the force and the cilium deformation in Fourier modes 
in time. For simplicity, we only consider here the first 
two Fourier components and do not look at the Fourier 
component h2 that characterizes the asymmetry of the 
beating. The Fourier components of the velocity are re- 
lated to the Fourier components of the force by 



E 



(fe,0#(o,o) 



(fc2 +?2)5/2 



(39) 



v (s ) ~ — 



Note that /C(0) — JC. The function JC{if) is plotted on 
Fig. IHIfor a lattice of 10^ cilia. 

Note that for two particular values of that we denote 
by and this function vanishes, IC(ips) = IC{(pa) = 0, 
as in the case where the relative phases of the cilia are 
randomly distributed. This corresponds to a constant 
flow with no oscillating component. 
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It is convenient to rewrite the external velocity as u{tf) 
iu)<Ci^{ip) with 



7M 



dshi{s)s (43) 



The constant Ci can be determined self-consistently as 
it varies linearly with hi. We obtain 



(44) 



with 



FIG. 8: /C((/p) over one period {ip G [0, 27r]). Some remark- 
able values: A^(0) = JCmax — 4.52; IC{tt) — tCmin — —2.32; 
/C((y3s) = IC{ipa) = with (y3s ~ 1.34 and (pa ~ 4.94. 



We now define the two dimensionless velocities U and 
u{<f) by 



U 



K 2Trrid^ 



3/C(y.)/,£ 



(40) 



The equations of motions of the Fourier components /iq 
and hi can then be written as 



/iQ — Kh^ = Us 
hi + Xihi + iiuhi = u{(p)s 



(41) 



In writing Eq. we only kept the term Us — 0{h^) 
that breaks the left-right symmetry and that lead to 
ho 7^ 0, ignoring any other term that would not create a 
macroscopic motion. 



2. Beating pattern and metachronal waves 

We first study the equation of motion of the first 
Fourier mode in Eq. 1411 which corresponds to the os- 
cillatory motion of the cilium. The right hand side of 
this equation of motion does not vanish due to the ex- 
istence of an oscillatory external flow due to the other 
cilia. Note however that we have not treated in details 
the hydrodynamic interactions for one cilium and that 
we have only taken them into account through the two 
local friction coefficients and ^|| . We are here more 
interested in the qualitative aspects of the coordination 
between cilia than in the accurate calculation of the flows 
created by each cilium. 

The general solution of Eq. ^]can be written as hi — 
h\ + hFi with 



h^{s) = Aie'i^' + Bie-'i^'' + Cie'i^' + Die-'^^' 

m = ^-s 



(42) 



l-7(^)/3 



1 



(45) 



The effect of the hydrodynamic interactions between cilia 
is embodied here in the coefficient 7(<^). The variation 
of this coefficient with the phase difference ip is similar 
to that of IC{(p). The limit where 7(<p) — 0, leads back to 
the previous situation were one cilium is beating alone; 
it may however correspond to the finite phase shifts be- 
tween cilia Lp — ips OT ipa- 

The four boundary conditions on hi can as before be 
written in a matrix form and the oscillation threshold and 
the beating frequency can be determined as the zeros of a 
determinant insuring the consistency of this matrix equa- 
tion. This leaves an unknown amplitude of the beating 
motion that could only be calculated by expanding the 
equation of motion to higher order. The beating pattern 
can then be written as 

hi{s) = Ai[£^{qi,s)+bi£^{-qi,s)+ci£^{q2,s)+di£^{-q2., 
with 



£v{q,s) 



= 9S 



■ Piq,(p)s 



The values of both the oscillation threshold ilc and the 
frequency Uc depend on the phase shift between cilia ip, 
through 7(93). We first discuss the variation of this bi- 
furcation point with the constant 7(1^9) which is a more 
convenient variable. On Fig. |3 we plot flc and the crit- 
ical frequency fc against 7. 

There is a local minimum of Qc for 7* ~ —1.15 and a 
local maximum for 7 ~ 1. The beating frequency fc, is 
a decreasing function of 7. 

We here need a selection criterion that determines 
the value of the phase shift between cilia. The sim- 
plest conjecture for the selection criterion is that the sys- 
tem chooses the local minimum of flc corresponding to 
7* ~ —1.15. This corresponds to a metachronal wave 
propagating in the assembly of cilia, as widely confirmed 
by experimental observations (^26. 27] for instance). 

With this selection criterion, the oscillation thresh- 
old is — 6.538 10~^ and the critical frequency is 
fc — 31Hz. The hydrodynamical couplings between 
cilia decrease the oscillation threshold flc and increase 
the critical frequency fc. The coordination between cilia 
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FIG. 9: Oscillation threshold flc and critical frequency fc as 
functions of 7 for 7 £ [—3,2]. flc has a local minimum that 
corresponds to the existence of metachronal waves. 




the parameters. If we take ij = iijw and d/L = 1 
so that our calculations remain consistent and in or- 
der to be close to what is observed experimentally, then 
IC{>f*) ~ -0.07 which yields ip* ~ ±1.37 ~ ±0.447r. This 
value corresponds to a wavelength A — A.6d ~ 5d for 
the metachronal waves or approximatively 6 cilia, which 
is the correct order of magnitude ( the wave length is 7 
cilia in ,27.]'). 



3. Global flow and left-right symmetry breaking 

We now discuss the left-right symmetry breaking and 
the appearance of a global flow. We solve Eq. ^]for the 
zeroth Fourier component of the deformation, with the 
same boundary conditions as before, in the limit K 0. 
We obtain 



ho(s) = U—(l-- + —) 



(46) 



The cilium oscillates around a curved average position 
/iQ 7^ if [7 ^ 0, if there exists a global flow. We show 
below that this is possible within a certain range of pa- 
rameters. 

We define the two dimensionless functions 



Hois) 



_ U 
hijs) 



6 



20^ 



Sipiqi, s) + bi£y,{~qi,s) 
Ci£cp{q2,s) + di£^(-q2, s) 



The determination of the average velocity U requires the 
calculation of the integral of /q defined in Eq. |23 we 
obtain 



(47) 



with 



ds '^[2HqHiHI - HI / ds'HoHi] (48) 



FIG. 10: Beating pattern of a cilium in an array in the pres- 
ence of a metachronal wave. The pattern is different from 
that of an isolated cilium mostly around the basis. The first 
Fourier component 25R[/iie'"*] at various time steps during a 
period is plotted. The parameters are Ai = 1/70; the maxi- 
mum deformation is hmax — 0.15. 



favors cilium beating by creating a metachronal wave cor- 
responding to 7 < 0. 

The beating pattern is slightly changed as shown on 
Fig. [TUl where we have plotted 23?[/iie'"*] at different 
time steps with the same amplitude Ai — 1/70. 

The phase difference ip* between two consecutive cilia 
corresponding to 7* ~ —1.15 depends on the values of 



which can be numerically calculated knowing Hq and Hi . 
dp depends on ip through Hi. Using the value of (j) cor- 
responding to metachronal waves, we obtain 

~ 34.5 

A self-consistent equation is then be obtained for the 
average velocity U 



U = 



■nrjd^ 



Aiuo U 



(49) 



If < this equation has the only solution U = i) and 
no global flow can exist, the left-right symmetry is not 
broken. If <C„ > this equation can have two extra non 
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zero solutions U ^ Q corresponding to a global flow along 
the X axis given by 

< V{s,t) K,(s) = Us 

and the left-right symmetry is then broken. 

The condition for appearance of a global flow is 

"TT > ^ 10"7 

TT rj K d'^ 

As for the oscillation amplitude, our calculation only 
gives the threshold of appearance of the global flow. A 
determination of the actual value of the velocity would 
require an expansion of the equations of motion to higher 
orders. 



V. DISCUSSION AND CONCLUDING 
REMARKS 

We have studied in this paper how hydrodynamic in- 
teractions between cilia contribute to the coordination 
of the beating motion in ciliated cells. Three major ef- 
fects have been studied, the spontaneous alignment of 
an array of cilia, the breaking of the symmetry of the 
beating and the appearance of a macroscopic flow and 
the existence of metachronal waves. We have shown for 
all these problems that there exist a dynamic transition 
where symmetry is broken and where a coordination be- 
tween the beating of neighboring cilia appears. 

Our work is based on several simplifying approxima- 
tions that we believe make the analysis tractable ana- 
lytically but that should preserve the essential physical 
effects. We only studied hydrodynamic interactions be- 
tween distant cilia that can be treated by introducing 
simple distribution of forces in the fluid to describe the 
motion of one ciiium. This is rarely true experimentally 
but the hydrodynamic interactions between closer cilia 
are even stronger and strongly favor the transitions that 
we study. We have replaced the complex architecture 
of the axoneme by two microtubules sliding against one 
another under the action of dynein motors which are de- 
scribed by a two state model for molecular motors as 
done earlier by Camalet and Jiilicher. This is a rather 
sketchy description but it allows a calculation of the in- 
ternal forces that drive the ciiium motion and it gives 
some physical insight. Future work will have to take into 
account the nine-fold symmetry of the axoneme and the 
influence of its central doublet. Finally, we have only 
considered small amplitude beating. This is sufhcient to 
determine the oscillation threshold but it does not allow 
a quantitative comparison between the calculated beat- 
ing and the experimental one that often occur far from 
any threshold. All our results are qualitatively consistent 
with the experimental observations and for example the 
beating frequency is close to both the experimental ones 
and to the ones obtained in numerical simulations [3^ . 

The essential result of our work is the natural emer- 
gence of metachronal waves and of a macroscopic flow 



created by an array of cilia if the amplitude of beat- 
ing is large enough. The criterion for appearance of the 
global uniform component of the flow given by Eq. I5UI 
requires only very small amplitudes {Ai > 5 10^^) which 
means that as long as the left-right symmetry is broken, a 
macroscopic flow should appear. An essential ingredient 
for the macroscopic flow to appear is that the constant 
Cip defined in Eq. ^Hlbe positive so that the average force 
created by one ciiium favors the flow and does not oppose 
it (which occurs ii (p = 0). 

As long as we allow a constant phase shift between 
neighboring cilia we observe metachronal coordination as 
a consequence of hydrodynamic interactions and of the 
internal beating mechanism of the ciiium. A selection 
criterion is then needed for these waves. We have con- 
jectured that the existing metachronal wave is the one 
that corresponds to the local minimum of the oscillation 
threshold. A more complete calculation that goes far be- 
yond the scope of this work would have to consider the 
nucleation of the metachronal wave and to determine the 
fastest growing wave. One of the interesting predictions 
of our calculation is that the existence of metachronal 
waves leads to a flow which is far more stationary than 
if all the cilia were beating in synchrony. The oscillat- 
ing component of the flow is proportional to the constant 
IC{lp) (see Eq. I4U|I which has a much smaller value when 
metachronal waves exist {IC{ip*) ~ —0.07) than if all cilia 
are beating in synchrony (/C(0) = /C ~ 4.52). Metachro- 
nism thus contributes to the creation of a very steady 
movement of swimming organisms that could for example 
make easier the detection of the organism environment. 

Our most important conclusion is the idea that 
metachronism and the existence of macroscopic flow 
around ciliated organisms can exist as self-organized phe- 
nomena driven by hydrodynamic couplings. We must 
stress however that other mechanisms could be at the 
origin of these cooperative effects. 
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Appendix I: Average force created by a single 
beating ciiium in a viscous fluid 

The aim of this appendix is to calculate the force and 
momentum averaged over one time period created by a 
general periodic beating of a single ciiium. We make 
two assumptions: the beating is planar and there is a 
stationary external flow. In section Hvl the average flow 
is created by the neighboring cilia. 

We call (f> the angle between the plane of beating and 
the direction of the external flow V that we take along 
the X axis. The ciiium of length L is located at the ori- 
gin and it is fixed at its basis. We denote by h{s,t) 
the distance between a point at arclength s on the cii- 
ium and the z axis at time t and by Z(s,t) the distance 
between a point at the arclength s on the ciiium and 
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FIG. 11: Sketch of a beating cilium in a plane at an angle 
with the direction of the external flow V. 



the xy plane. The angle between the tangent vector t 
to the cilium and the z axis is denoted by ipls^t) (see 
Fig. Ill() . The coordinates of the tangent vector are 
t = (cos (j) sin ip, sin (f) sin ip^ cos ip). The angle "0 is related 
to the cilium deformation h by simjj = dgh 

The point on the cilium at the arclcngth s is located 
at position X = {x,y, z), with: 



X = cos(j) / du simp (u, t) ~ h{s,t) cos (j) 
Jo 

y = simp / du simp{u,t) = h{s,t) sine/) 
Jo 



du cos'0(u,t) = Z{s,t) 



The velocity of this point is calculated by derivation with 
respect to time, v = dtX. 

The force per unit length exerted by the cilium on the 
fluid expressed in the Frenet basis (t, n, b) is: 



/ = {^\\tt + Un n + Ub h){v- V) 



(51) 



where ^|| and are the two local friction coefficients for 
tangential and normal motion respectively. We decom- 
pose this force as a sum of two forces, depending 
on the local velocity and //'"^ depending on the exter- 
nal flow velocity and calculate the average force over a 

1 T ^ 
beating period < f > ^ 7f L dt f(t). 



calculated 






< fl^"' > = 


(a- 


. cos (p 


^ cheat ^ 




sin0 
2 


< > = 






A(?/,s) = 




- V'(^^) 



du < dtip{u) cos A(u, s) > 
du < dttp{u) cos A(w, s) > 



(52) 



This force is proportional to (^j_ — ^||) as mentioned in 
sectional The difference between the two local friction 
coefficients ^± and ^|| is at the basis of the flow generation 
by an assembly of beating cilia. Indeed, this is because 
the shape of the beating in the effective stroke is differ- 
ent from that in the recovery stroke that a force can be 
exerted in the fluid on average. 

The average force due to the external flow is 



<fy 



flow ^ _ 
flow ^ _ 



(C± - £.\\)Vcos^ (P < sin^ ip> -£.±V 
{£,± — ^||)V^cos0sin(/) < sin^ tp > 
(C_L — < sini/)C0Si/' > 



The average beating force < f^'^'^* > can be explicitly 



It is important to note that < fl^°^{s) >< : this force 
opposes the flow. The last term of //'"^ is a static term, 
whereas the first positive term depends on the beating 
pattern and reduces the effects of this static term. In 
an assembly of cilia, the external velocity is due to the 
beatings of the other cilia which are themselves created 
by the forces on these cilia. 

In sectional we introduce a viscous coefficient a which 
characterizes the tendency for a cilium, beating in a 
plane at an angle cp with the flow, to align with the 
other cilia. A torque along the z axis due to the flow 

j^jflow ^ _^JJ 

sin <p is exerted on this cilium. We now 
express a as a function of the cilium beating pattern. We 
call ruz = — (X X f).ez the torque along z exerted by the 
fluid on the cilium per unit of length (the minus sign is 
due to the fact that / is the force exerted by the cilium 
on the fluid). The local torque per unit length exerted 
by the fluid on the cilium is 

mz{s, t) — —^±Vh(s, t) sin p = —^±Lh{s, t)Z{s, t)U sin (p 

where we have used the dimensionless coordinates s = 
s/L, h = h/L. The total momentum along z averaged 
over time, is obtained by integration 

Mz = -£,±L^ [ ds <h{s,t)Z{s,t) > UsintP (53) 
Jo 

This defines the friction coefficient a: 

a^^^L^ [ ds <h{s,t)Z{s,t) > (54) 

which can be calculated if the motion of the cilium is 
known. 
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Appendix II 

In this appendix, we derive the equations satisfied by 
the Fourier components of the deformation /i of a single 
beating ciHum and we determine the threshold of spon- 
taneous oscillations of the cilium. 



Our choice of a symmetric potential Wi imposes that 
a change A ^ — A must change / ^ — /. This sym- 
metry imposes thus /(^'^') — 0. The only non-vanishing 
coefficient at linear order is f^^^ — nuj)Sn.i with 



x(r2, Lu) — — Xiuj + 



-^57) 



Fourier mode expansion 

Axoneme beating is periodic and can be studied by 
expansion in Fourier modes in time of all the physical 
parameters: 



h{s,t) ^ ^ hn{s)e 



n— — oo 



The definition is similar for the other parameters. Start- 
ing from Eq. 1171 and Eq. 1141 we obtain the Fourier com- 
ponents: 

' Jo 



Ao. 



1 



k + inujj Jq 



dsf ni-s) 



(55) 



In order to determine the non linear relationship be- 
tween / and A, we follow the lines of [s^ and write: 

fn - fi'^ + E /IP A. + E /(^A^A™ + O(A^) 



Im 



The force and the sHding displacement are thus related 
by 

U = xin,nLu)A^+0{\A\^) ^ x{^,riLu){Aon+ak,,+0{\h\^)) 

(58) 

This relationshipl58lmodels the response of the molecular 
motors to the bending of the axoneme. 
From Eq. 1551 we obtain 



Aon — 



x(r2, nuj)a 



k + inuj — x(^; nuj)L 



K{L)+0{\h\^) (59) 



The coefficients fi%,...,nk can be calculated by first Fourier component reads 
rewriting Eq. 1151 as 



We solve the equation of motion of the cilium fEa ll8|l for 
each order of the Fourier expansion. 



Equation of motion of the Fourier modes 

We look for an approximate solution of the form 

h{s, t) ~ ho{s) + /ii(s)e*'^' + h2{s)e^''^' + c.c. 

At linear order, there is no coupling between the modes 
and using Ea ll8l the equation of motion of the n"" 



P„ = R5n,0 - ^\ E l'Sn,l+mAld^Pm (56) 

Ini 



where 



UJ2{0 l + f}sin2«/0 

R = p ■ = p- 



1 + 



is the static probability {lu — 0), corresponding to a 
medium with not enough ATP to generate the beating. 
Inserting the ansatz 

Pn = RSnfl + E pIu^I + E pIIAi^-^ + 0{A^) 
I Irn 

into Eq. 1561 we obtain a recursion relation for the 

-i n,ni,...,nfc . 



.... x{^,nuj)a'^ ■■ .nuj£,j_ 
hn-\ hn + I K = (60) 

K K 

It is convenient to introduce the dimensionless variables: 



- /r - - -fa m ^ 

S = S/ L LU = LU Xn = X(": '^'^j = '^'-^j 

K K 

In dimensionless form, Eq. 1571 can be written 



x{n,Lu) = ~K -Xiuu+'^pU * 



2 V + iLU 



with 



2t2 



K 



■K A 



A D 



V U ■ 



■U 



that now allows us to calculate /n*li,...,n 



We have anticipated here the fact that il <C 1. 

Defining, h — h/L, and denoting by a dot the deriva- 
tion with respect to s, we obtain the equation of motion 
Eq. 12 II and the boundary conditions given by Eq. I22I23I 
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Beating motion 



with 



In the absence of external flow only the first Fourier 
component of h does not vanish and satisfies the equation 
of motion : 

h 1 + xi hi + iuj hi = (61) 
where the relevant dimensionless parameters are 



J'{q)=ei{q'+xiq + ri) 



Xi = xi^,^) 



Ti = 



xi 



k — Xi + ^1^ 



The boundary conditions are given by Eq. 1^ and ESI for 
n = 1. The solution to this linear equation is a superposi- 
tion of exponentials given by Eq. 1241 The four boundary 
conditions on hi can be written in a matrix form: 



Mi(f7,6j).Ai = 



(62) 



The systemic has non trivial solutions only if 



detMi(f7,a;) = 



(63) 



where Ai is the vector made by the amplitudes of the 
exponentials in Eq. 1241 and the matrix Mi is given by 



Mi(fi,u;) 



1 



1 



1 



1 



qi -qi q2 -92 

T{qi) T{-qi) T{q2) T{-q2) 



qfe^^ 



qje-_^ 



Since Eq. 1631 is a complex equation, it determines both 
the oscillation threshold and the dimensionless beat- 
ing frequency uic- 



[1] B. Alberts, D. Bray, J. Lewis, M. Raff, K. Roberts, and 

J. Watson, Molecular Biology of the Cell (Garland, New 

York, 1994). 
[2] R. Anderson, J. Cell. Biol. 54, 246 (1972). 
[3] M. Sleigh, The Biology of Cilia and Flagella (Pergamon 

Press, Oxford, 1962). 
[4] H. Hagiwara, N. Ohwada, and K. Takata, Int. Rev. Cytol. 

234, 101 (2004). 
[5] M. Porter and W. Sale, J. Cell. Biol. 151, F37 (2000). 
[6] M. Wargo and E. Smith, Proc. Natl. Acad. Sci. USA. 

100, 137 (2003). 
[7] S. Camalet, Ph.D. thesis, Universite Paris VI, Paris 

(2001). 

[8] S. Camalet and F. Jiilicher, New J. Phys. 2, 24.1 (2000). 

[9] Y. Naitoh and H. Kaneko, Science. 176, 523 (1972). 
[10] Y. Naitoh and H. Kaneko, J. Exp. Biol. 58, 657 (1973). 
[11] I. Gibbons, J. Biophys. Biochem. Cytol. 11, 179 (1961). 
[12] S. Sorokin, J. Cell. Sci. 3, 207 (1968). 
[13] E. Boisvieux-Ulrich, M. Laine, and D. Sandoz, Biol. Cell. 

55, 147 (1985). 
[14] B. Afzelius, Int. Rev. Exp. Path. 19, 1 (1979). 
[15] N. Liron and S. Mochon, J. Fluid. Mech. 75, 593 (1976). 
[16] J. Blake, Proc. Camb. Phil. Soc. 70, 303 (1971). 
[17] J. Beisson and M. Jerka-Dziadosz, Biol. Cell. 91, 367 
(1999). 

[18] I. Gradshteyn and I. Ryzhik, Table of integrals, series 
and products (Academic Press, Orlando, 1980). 

[19] J. M. Borwein and P. B. Borwein, Pi and the ACM: 
A Study in Analytic Number Theory and Compu- 



tational Complexity (Wiley, New York, 1987), 291. 
http: / /mathworld.wolfram.com/DoubleSeries.html. 
[20] A. Hilfinger, Ph.D. thesis, Technische Universitat, Dres- 
den (2005). 

[21] G. Vernon and D. WooUey, Biophys. J. 87, 3934 (2004). 
[22] J. Prost, J.-F. Chauwin, L. Peliti, and A. Adjari, Phys. 

Rev. Lett. 72, 2652 (1994). 
[23] F. Jiilicher and J. Prost, Phys. Rev. Lett. 75, 2618 

(1995). 

[24] S. Ishijima and Y. Hiramoto, Cell Struct. Funct. 19, 349 
(1994). 

[25] M. Murase, The dynamics of cellular motility (Wiley, 

Chichester, 1992). 
[26] E. Horstmann, Movie, IWF, Gottingen (1959), 

http://mkat.iwf.de. 
[27] H. Machemer, J. Exp. Biol. 57, 239 (1972). 
[28] S. Gueron and K. Levit-Gurevich, Biophys. J. 74, 1658 

(1998). 

[29] P. Sugrue, M. Hirons, J. Adam, and H. M.E., Biol. Cell. 
63, 127 (1988). 

[30] M. Hayashi, T. Yagi, Y. K., and R. Kamiya, CeU. Motil. 

Cytoskeleton. 41, 49 (1998). 
[31] M. Cosentino Lagomarsino, B. Bassetti, and P. Jona, 

Eur. Phys. J. B. 26, 81 (2002). 
[32] S. Gueron, K. Levit-Gurevich, N. Liron, and J. Blum, 

Proc. Natl. Acad. Sci. USA. 94, 6001 (1997). 
[33] F. Julicher and J. Prost, Phys. Rev. Lett. 78, 4512 

(1997). 



